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ENTROPY-DISSIPATING SEMI-DISCRETE RUNGE-KUTTA SCHEMES 
EOR NONLINEAR DIFFUSION EQUATIONS 


ANSGAR JUNGEL AND STEFAN SCHUCHNIGG 

Abstract. Semi-discrete Runge-Kutta schemes for nonlinear diffusion equations of para¬ 
bolic type are analyzed. Gonditions are determined under which the schemes dissipate the 
discrete entropy locally. The dissipation property is a consequence of the concavity of the 
difference of the entropies at two consecutive time steps. The concavity property is shown 
to be related to the Bakry-Emery approach and the geodesic convexity of the entropy. The 
abstract conditions are verified for quasilinear parabolic equations (including the porous- 
medium equation), a linear diffusion system, and the fourth-order quantum diffusion equa¬ 
tion. Numerical experiments for various Runge-Kutta finite-difference discretizations of 
the one-dimensional porous-medium equation show that the entropy-dissipation property 
is in fact global. 


1. Introduction 

Evolution equations often contain some structural information reflecting inherent phys¬ 
ical properties such as positivity of solutions, conservation laws, and entropy dissipation. 
Numerical schemes should be designed in such a way that these structural features are pre¬ 
served on the discrete level in order to obtain accurate and stable algorithms. In the last 
decades, concepts of structure-preserving schemes, geometric integration, and compatible 
discretization have been developed [7], but much less is known about the preservation of 
entropy dissipation and large-time asymptotics. 

Entropy-stable schemes were derived by Tadmor already in the 1980s [20] in the context 
of conservation laws, thus without (physical) diffusion. Later, entropy-dissipative schemes 
were developed for (hnite-volume) discretizations of diffusion equations in [21 HOI El] • In [5] , 
a hnite-volume scheme which preserves the gradient-how structure and hence the entropy 
structure is proposed. All these schemes are based on the implicit Euler method and are of 
hrst order (in time) only. Higher-order time schemes with entropy-dissipating properties 
are investigated in very few papers. A second-order predictor-corrector approximation was 
suggested in [19], while higher-order semi-implicit Runge-Kutta (DIRK) methods, together 
with a spatial fourth-order central hnite-diherence discretization, were investigated in [5]. 
In [H EZ], multistep time approximations were employed, but they can be at most of 
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second order and they dissipate only one entropy and not all functionals dissipated by the 
continuous equation. In this paper, we remove these restrictions by investigating time- 
discrete Runge-Kutta schemes of order p > 1 for general diffusion equations. 

We stress the fact that we are interested in the analysis of entropy-dissipating schemes 
by “translating” properties for the continuous equation to the semi-discrete level, i.e., we 
study the stability of the schemes. However, we will not investigate convergence, stiffness, 
or computational issues here (see e.g. m)- 

More precisely, we consider time discretizations of the abstract Cauchy problem 

(1) dtu{t) + A[u{t)] =0, t > 0, m(0) = 

where A : D{A) —)■ X' is a (differential) operator dehned on D{A) C X and X is a 
Banach space with dual X'. In this paper, we restrict ourselves to diffusion operators 
A[u] dehned on some Sobolev space with solutions u : x (0, cxd) —)■ M”, which may be 
vector-valued. A typical example is A[u\ = div(a(M)VM) dehned on X = L^(fl) with 
domain D{A) = where a : M —?■ M is a smooth function (see section [3]). Equation 

([I]) often possesses a Lyapunov functional H[u\ = J^h{u)dx (here called entropy), where 
h : M"' —)■ M, such that 

[m] = f h'{u)dtudx = — f h'{u)A[u]dx < 0, 
dt Jq Jn 

at least when the entropy production J^h'{u)A[u]dx is nonnegative. Here, h' is the de¬ 
rivative of h and h'{u)A[u] is interpreted as the inner product of h'{u) and A[u] in M”'. 
Furthermore, if h is convex, the convex Sobolev inequality j^h'{u)A[u]dx > kH[u] for 
some A > 0 may hold [6], which implies that dH/dt < —kH and hence exponential con¬ 
vergence of H[u] to zero with rate k. The aim is to design a higher-order time-discrete 
scheme which preserves this entropy-dissipation property. 

To this end, we propose the following semi-discrete Runge-Kutta approximation of (fT|): 
Given G X, dehne 


( 2 ) 


i=l 


-A 


u 


k-l 


i=i 


CLijKj 


i = l,...,s. 


where t^ are the time steps, r = t^ — t^~^ > 0 is the uniform time step size, approximates 
u{tA), and s > 1 denotes the number of Runge-Kutta stages. Since the Cauchy problem is 
autonomous, the knots Ci,... ,Cs are not needed here. In concrete examples (see below), 
are functions from to M”. If aij = 0 for j > i, the Runge-Kutta scheme is explicit, 
otherwise it is implicit and a nonlinear system of size s has to be solved to compute Ki. 
We assume that scheme ([2]) is solvable for : Q ^ M"". 

Given h : M"" —)■ M, we wish to determine conditions under which the functional 

Hlu’^] = f h{u^{x))dx 

Jn 


( 3 ) 
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is dissipated by the numerical scheme ([2]), 

(4) H[u’^] +T [ A[u’^]h'{u^)dx < H[u’^-% keN. 

Jn 

In many examples (see below), f^A[u’‘]h'{u^)dx > 0 and thus, the function k i-A H[u^] is 
decreasing. Such a property is the hrst step in proving the preservation of the large-time 
asymptotics of the numerical scheme (see Remark [2]). 

Our main results, stated on an informal level, are as follows: 

(i) We determine an abstract condition under which the discrete entropy-dissipation 
inequality (jl]) holds for sufficiently small > 0. This condition is made explicit for 
special choices of A and h, yielding entropy-dissipative implicit or explicit Runge- 
Kutta schemes of any order. 

(ii) Numerical experiments for the porous-medium equation indicate that may be 
chosen independent of the time step k, thus yielding discrete entropy dissipation for 
all discrete times. 

(iii) We show that for Runge-Kutta schemes of order p > 2, the abstract condition in (i) 
is exactly the criterion of Liero and Mielke [18] to conclude geodesic 0-convexity of 
the entropy. In particular, it is related to the Bakry-Emery condition [T]. 

Let us describe the main results in more detail. We recall that the Runge-Kutta scheme 
d^j) is consistent if ^^^1 X]i=i h = ^- Furthermore, if X]i=i = |, h is at 

least of order two [HI Chap. II]. We introduce the number 

S 

(5) Crk = 2 bi{l — Cj), 

i=l 

which takes only three values: 

Crk = 0 for the implicit Euler scheme, 

Crk = 1 for any Runge-Kutta scheme of order p >2, 

Crk = 2 for the explicit Euler scheme. 

The first main result is an abstract entropy-dissipation property of scheme ([2D for en¬ 
tropies of type (jsp. 

Theorem 1 (Entropy-dissipation structure I). Let h G let A : D{A) —)■ X' he 

Frechet differentiable with Frechet derivative DA[u] : X ^ X' at u E D{A), and let (pA) 
he the Runge-Kutta solution to ([2D. Suppose that 

(6) := [ {CnKh'{u’^)DA[u’^]{A[u'^]) + h"{u’^){A[u’^]fi)dx > 0. 

Jn 

Then there exists > 0 such that for all 0 < r < , 

(7) H[u’^] +T [ A[u’^]h\u'^)dx < H[u^-^]. 

Jn 
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We assume that the solutions to ([2]) are sufficiently regular such that the integral ([6]) 
can be dehned. In the vector-valued case, h'\u^) is the Hessian matrix and we interpret 
h”{u^){A[u^]Y as the product A[u^Yh”{u^)A[u^]. For Runge-Kutta schemes of order p > 
2 (for which Crk = 1), fhe integral dH]) corresponds exactly to the second-order time 
derivative of H[u{t)] for solutions u{t) to the continuous equation ([1]). Observe that the 
entropy-dissipation estimate (jTj) is only local, since the time step restriction depends on the 
time step k. For implicit Euler schemes (and convex entropies h), it is known that can be 
chosen independent of k. For general Runge-Kutta methods, we cannot prove rigorously 
that stays bounded from below as k ^ oo. However, our numerical experiments in 
section [3 indicate that inequality ([7]) holds for sufficiently small r > 0 uniformly in k. 

Remark 2 (Exponential decay of the discrete entropy). If the convex Sobolev inequality 
A[u’‘]h'{u^)dx > hiH[u^] holds for some constant a > 0 and if there exists r* > 0 such 
that > T* > 0 for all fc G N, we infer from ([7]) that for r := r*, 

< (1 + = exp{—riKt’^)H[u°], where p = 

KT 

which implies exponential decay of the discrete entropy with rate pn. This rate converges 
to the continuous rate a as r —)■ 0 and therefore, it is asymptotically sharp. □ 

Theorem [Dean be generalized to a larger class of entropies, namely to so-called first-order 
entropies 

(8) F[u'^]= [ \Vf{u^)\^dx, 

Jn 

where, for simplicity, we consider only the scalar case with / ; M —)■ M. An important 
example is the Fisher information with f{u) = y/u. 

Theorem 3 (Entropy-dissipating structure H). Let f E C^(M), let A : D{A) -E X' be 
Frechet differentiable, and let {vf) be the Runge-Kutta solution to (j2]). Assume that the 
boundary condition X f{u^) ■ u = 0 on dfl is satisfied. Furthermore, suppose that 

Ii--= [ (|V(/'(n^)A[n^]|2-CRKA/(n^)/'(n^)DA[n'=](A[n^]) 

(9) ^ 

Then there exists > 0 such that for all 0 < r < r^, 

F[u'‘]+t f A[u^]f{u^)dx < F[u’^-^]. 

Jn 

The key idea of the proof of Theorem [1] (and similarly for Theorem [3]) is a concavity 
property of the difference of the entropies at two consecutive time steps with respect to 
the time step size r. To explain this idea, let u := be hxed and introduce v{t) := u^~^. 
Clearly, n(0) = u. A formal Taylor expansion of G(r) := H[u] — H[v{t)] yields 

H[u’^] - H[u^-^] = G{t) = G(0) + rG'(O) + yG"(e^), 



ENTROPY-DISSIPATING SEMI-DISCRETE RUNGE-KUTTA SGHEMES 


5 


where 0 < .^^ < r. A computation, made explicit in section [2l shows that G'(0) = 
J^A[u’^]h'{u’^)dx and G"(0) = —Jq. Now, if G"(0) < 0, there exists > 0 such that 
G"(r) < 0 for r e [0,r^] and in particular G"(^^) < 0. Consequently, G(r) < rG'(O), 
which equals (jlj). The dehnition of u(r) assumes implicitly that ([2]) is backward solvable. 
We prove in Proposition [5] below that this property holds if the operator A is a smooth 
self-mapping on X. 

Remark 4 (Discussion of r^). Since {u^) is expected to converge to the stationary solution, 
lim^^oo Iq = 0. Thus, in principle, for larger values of /c, we expect that becomes 
smaller and smaller, thus restricting the choice of time step sizes r. However, practically, 
the situation is better. For instance, for the implicit Euler scheme, if h is convex, we obtain 

H[u^] - H[u’^-^] < f -u^-^)dx = -r f h’{u^)A[u^]dx 

Jn Jn 

for any value of r > 0. Moreover, for other (higher-order) Runge-Kutta schemes, the 
numerical experiments in section [7| indicate that there exists r* > 0 such that G"(r) < 0 
holds for all r G [0, r*] uniformly in /c G N. In this situation, inequality ([7]) holds for all 
0 < r < r*, and thus our estimate is global. In fact, the function G" is numerically even 
nonincreasing in some interval [0, r*] but we are not able to prove this analytically. □ 

The second main result is the specihcation of the abstract conditions ([6]) and ([9]) for 
a number of examples: a quasilinear diffusion equation, porous-medium or fast-diffusion 
equations, a linear diffusion system, and the fourth-order Derrida-Lebowitz-Speer-Spohn 
equation (see sections [MSI for details). For instance, for the porous-medium equation 

dtu = A(m^) in D, t > 0, • z/ = 0 on dQ, ■u(O) = 

we show that the Runge-Kutta scheme scheme satishes 

+t( 3 f {u^Y+^-^\Vu^\^dx < H[u^-\ where H[u] = —^ 

Jn a{a + ^)Jn 

for 0 < r < and all {a, (3) belonging to some region in [0, oo)^ (see Figure [T] below). For 
a = 0, we write H[u] = J^u(\ogu — l)dx. In one space dimension and for Runge-Kutta 
schemes of order p>2, this region becomes —2 < a — (3 < 1, which is the same condition 
as for the continuous equation (except the boundary values). Furthermore, the hrst-order 
entropy (IH]) is dissipated for Runge-Kutta schemes of order p > 2, in one space dimension, 

F[u’^] + f < F[u’^-^], where F[u] = [ 

Jn Jn 

for 0 < r < and all (a, (3) belonging to the region shown in Figure [2] below, and Ca,p > 0 
is some constant. This region is smaller than the region of admissible values {a, 13) for the 
continuous entropy. The borders of that region are indicated in the hgure by dashed lines. 

The proof of the above results, and namely of G"(0) < 0, is based on systematic in¬ 
tegration by parts [H]. The idea of the method is to formulate integration by parts as 
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manipulations with polynomials and to conclude the inequality G"(0) < 0 from a polyno¬ 
mial decision problem. This problem can be solved directly or by using computer algebra 
software. 

Our third main result is the relation to geodesic 0-convexity of the entropy and the 
Bakry-Emery approach when Grk = 1 (Runge-Kutta scheme of order p > 2). Liero and 
Mielke formulate in [18] the abstract Cauchy problem ([1]) as the gradient flow 

dtu = —K[u]DH[u], t > 0, m(0) = 


where the Onsager operator K[u] describes the sum of diffusion and reaction terms. For 
instance, if A[u] = div(a(M)VM), we can write A[u] = dw{a{u)h''{u)~^Vh'{u)) and thus, 
identifying h'{u) and DH[u], we have K[u\^ = div(a(M)h"(M)“^V^). It is shown in [T8| 
that the entropy H is geodesic A-convex if the inequality 

( 10 ) := 

holds for all suitable u and We will prove in section |2] that 

G"(0) = 


Hence, if G"(0) < 0 then flTU]) with A = 0 is satished for u = and ^ = h'{u^), yielding 
geodesic 0-convexity for the semi-discrete entropy. Moreover, if G"(0) < — AG'(O) then we 
obtain geodesic A-convexity. Since G'(0) = —dH[u]/dt and G"(0) = —d‘^H[u]/dt‘^ in the 
continuous setting, the inequality G"(0) < —AG'(O) can be written as 


d^H 

df^ 


\u\ 




which corresponds to a variant of the Bakry-Emery condition [T], yielding exponential 
convergence of H[u] (if > r* > 0 for all k). Thus, our results constitute a hrst step 
towards a discrete Bakry-Emery approach. 

The paper is organized as follows. The abstract method, i.e. the proof of backward 
solvability and of Theorems [H and [3l is presented in section [2l The method is applied 
in the subsequent sections to a scalar diffusion equation (section |3|), the porous-medium 
equation (section |4|), a linear diffusion system (section |5|), and the fourth-order Derrida- 
Lebowitz-Speer-Spohn equation (section |6|). Finally, section [7] is devoted to some numerical 
experiments showing that G" is negative in some interval [0,r*]. 


2. The abstract method 

In this section, we show that the Runge-Kutta scheme is backward solvable if H is a 
self-mapping and we prove Theorems [D and El 

Proposition 5 (Backward solvability). Let G [0, cxo) x X, where X is some Banach 

space, and let A G G^(X, X) he a self-mapping. Then there exists tq > 0, a neighborhood 
V <Z X of u^, and a function v G G^([0,ro);X) such that ([2|) holds for u’‘~^ := v{t). 
Moreover, 

( 11 ) 


n(0) = 0, T^O) = A[v\, and n^^(O) = C^^kDA[u]{A[u]). 
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The self-mapping assumption is strong for differential operators A but it is somehow 
natural in the context of Runge-Kutta methods and valid for smooth solutions. 

Proof. The idea of the proof is to apply the implicit function theorem in Banach spaces (see 
[H Corollary 15.1]). To this end, we set u := and dehne the mapping J = (Jq, ..., J^) : 
M X ^ by 

y) = v-u + T'^ biki, where y = (fci, ...,ks,v), 


2=1 


Ji{T,y) = ki + A 


v + r'^ttijkj 

j=i 




The Frechet derivative of J in the direction of (r/i, yt), where yh = {khi, ■ ■ ■, khs, Vh), reads 
as 

s s 

DJo{T,y){Th,yh) =Vh + Th '^biki + r bi kfii , 

i=\ 

DJiir, y){Th, Vh) = hi + DA 


2=1 


v + T^ ttijkj ivh + Th'^ aijkj + D dijkhj j , 

L j = l J V j = l j = l / 

where f = 1,..., s. Let Tq = 0 and yo = {—A[u ],..., —A[u],u). Then J(ro, yo) = 0 and 

DJQ{To,yo){Q,yh) =Vh, DJi{To,yo){Q,yh) = kih + DA[u]{vh), i = l,...,s. 

The mapping yh i-A HJ(ro, |/o)(0, |/h) is clearly an isomorphism from onto By 

the implicit function theorem, there exist an interval U C [0, tq), a neighborhood V C 
of yo, and a function (k,v) G C^([0,ro); V) such that (fc,n)(0) = {—A[u ],..., —A[u],u) and 
J(r, /c(r), n(r)) = 0 for all r G [0, Tq). 

Implicit differentiation of J(r, fc(r), n(r)) = 0 yields 


0 = v'{t) bikiir) + D hKir), 


2=1 


2=1 


0 = k'AT) +DA 


i=i 


+ r ^ Oij /cj (r) ( n'(r) + ^ aij /cj (r) + r ^ /c' (r) ), 


i=i 


i=i 


where i = 1,..., s and r G [0, tq). Using h = ^ and = c*, we infer that 

S S 

n'(0) = -y^6i/ci(0) = '^biA[u] = A[u], 


2=1 


2=1 


K{^) — —DA[u] I A[u] — '‘^aijA[u] j = —(1 — Ci)DA[u]{A[u]). 

i=i 


( 12 ) 
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Differentiating Jo(r, /c(r), n(r)) = 0 twice leads to 


0 = n"(r) + 2 ^ biklir) + t ^ &i/c"(r). 

2=1 2=1 

Because of flT^ . this reads at r = 0 as 

S S 

v'\0) = -2^6i/c'(0) = 2'^bi{l - Ci)DA[u]{A[u]) = CkkDA[u]{A[u\). 


2=1 


2=1 


This hnishes the proof. 

We prove now Theorems [1] and [3l 


□ 


Proof of Theorem [3 We set u := u^. By Proposition [5], there exists a backward solution 
V G C^([0,ro)) such that r(0) = m, r'(0) = A[u], and r"(0) = C^kD 2\.[u]{A[u]). Further¬ 
more, the function G{t) = f^(h(u) — h(v(T)))dx satishes G(0) = 0, 

G'(0) = — f h'{v{0))v'{0)dx = — [ h'{u)A[u]dx, 

Jn Jn 

G"(0) = - f (h'(r(0))r"(0) + h"{v{0))v'{0y)dx 
Jn 

= — f {u)DA[u\{A[u\) + h"{u){A[u\)‘^)dx = —Iq < 0, 

Jn 

using the assumption. By continuity, there exists 0 < < Tq such that G'Jf) < 0 for 

0 < .^ < r^. Then the Taylor expansion G(r) = G(0) -|- G'(0)r -|- 

concludes the proof. □ 

Proof of Theorem 0 Following the lines of the previous proof, it is sufficient to compute 
G'(0) and G"{0), where now G{t) = /^(|V/(m)P — \ Vf{v{T))\‘^)dx. Using integration by 
parts and the boundary condition V/(r) ■ z/ = 0 on dfl, we compute 

G'(0) = - [ V/(r(0))-V(f(r(0)K(0))dx= [ Af{u)fiviT))A[u]dx, 

Jn Jn 

since r(0) = u and r'(0) = A[u]. Furthermore, again integrating by parts, 

G"(r) = - (|V(f(r(r))U(r))|V V/(r(r)) ■ V(/"(r(r))(r'(r))2) 

+ V/(r(r)) ■ V(/'(r(r))r"(r)) jdx 

I V(f (r(r))U(r)) 1^ - A/(r(r))/"(r(r))(U(r))2 


J^f{v{T))f{v{T))v"{T)) dx. 
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Since f"(0) = CYiKDA[u]{A[u]), this reduces at r = 0 to 

G"(0) = - f (|V(/'(«)>1H)P - Af(u)f"(u)(Alu])'^ - CnK^^f(u)f(u)DAlu](Alu])yx. 

This expression equals —Ii, and the result follows. □ 

Finally, we show that G"(0) for entropies (jS]) is related to the geodesic convexity condition 

of dg. 

Lemma 6. Let A[u] = K{u)DH[u] for some symmetric operator K ; D{A) —)■ X and 
Frechet derivative DH[u], let G he defined as in the proof of TheoremUl for a solution 
to the Runge-Kutta scheme ([2]) of order p > 2, and let M{u,f) be given by flTU]) . Then 

G"(0) = -2M{u’^,DH[u^]). 

Proof. The proof is just a (formal) calculation. Recall that for Runge-Kutta schemes of 
order p > 2, we have Crk = 1- Set u := and identify DH[u] with ^ = h'{u). Inserting 
the expression DA[u]{v) = DK[u]{y)h'{u) -|- K[u\h''{u)v into the dehnition of G"(0), we 
hnd that 

-G"(0) = {f,DA[u]{A[u])) + {A[u],h%u)A[u]) 

= DK[u]{A[u])f + K[u]h''{u)A[u]) + {A[uf h''{u)A[u]) 

= {f,DK[u]{K[u]i)f) + {f,K[u]h'\u)K[u]f) + {K[u%h" {u)K[u]i) 

= (e, DK[u]{K[u]m + 2(e, K[u]h"{u)K[u]f), 
since K[u] is assumed to be symmetric. Rearranging the terms, we obtain 

-G"(0) = + - {i,DK[u]{K[u]0) 

= 2{^,DA[u]iK[u]00 - {^,DK[u]iA[u])) = 2M(w,0, 

which proves the claim. □ 


3. Scalar diffusion equation 

In this section, we analyze time-discrete Runge-Kutta schemes of the diffusion equation 
(13) dtu = div(a(M)VM), f > 0, m(0) = 

with periodic or homogeneous Neumann boundary conditions. This equation, also includ¬ 
ing a drift term, was analyzed in [18] in the context of geodesic convexity. Our results are 
similar to those in [T8| but we consider the time-discrete and not the continuous equation 
and we employ systematic integration by parts [IT] . 

Setting p ( m ) = a{u)/h"(u) , we can write the diffusion equation as a formal gradient flow: 

dtU = —:= div{p.{u)'Vh'{u)), t > 0. 

We prove that the Runge-Kutta scheme ([2]) dissipates all convex entropies subject to some 
conditions on the functions /i and h. 
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Theorem 7. Let Q G be convex with smooth boundary. Let (u^) be a sequence of 
(smooth) solutions to the Runge-Kutta scheme ([2]) of the diffusion equation ffT5D . Let 
A; G N 6e fixed and be not equal to the constant steady state of (fT3l) . We suppose that 
for all admissible u, it holds that a{u) > 0, h"{u) > 0, 

(14) 6(m) := ^(Crk + 1) / fi{v)fi\v)h”{v)dv >0, 

^ JuQ 

(15) ^K«) < (C^RK + lW{u)fi{u)f 

d 

(16) (C'rk + 2 ) p,{u)p,"(u) + (Crk — l)p'{uy < 0. 

Then there exists > 0 such that for all 0 < t < , 

+ t [ h''{u^)aiu'^)\Vu^\‘^dx < H[u^-^]. 

Jn 

Conditions correspond to (4.12) in [18]. Condition (1T6|) is satisfied for concave 

functions p, except for the explicit Euler scheme (Crk = 2) for which we need additionally 
App" + {p'Y < 0. For the implicit Euler scheme, we may allow even for nonconcave 
mobilities p, e.g. p{u) = u'’' for 1 < 7 < 2 . 

Proof. According to Theorem [H we only need to show that Jq = —G"(0) > 0. To simplify, 
we set u := . First, we observe that the boundary condition Vu • = 0 on implies 

that 0 = dtVu-v = VdtU-v = —VA[m] -v on dVt. Using DA[u]{A[u]) = div{a'{u)A[u]Vu + 
a{u)VA[u]) = A{a{u)A[u]), the abbreviation = h'{u), and integration by parts, we 
compute 

G"(0) = — j {c-^iKh'{u)A{a{u)A[u]) + h'\u)(^di'v{p{u)Vh'{u))Y^dx 

CRKVh'(M) ■ V(a(M)A[M]) - h'\u){p\u)Vu ■ Vh\u) + p{u)Ah\u)Y^dx 

= -J (^CRKAfa{u)A[u] + h"{u) + p{u)A^^ "^dx. 

The boundary integrals vanish since Vu • z/ = VA[m] • = 0 on dQ. Replacing A[u] by 

div{p{u)'Vf) = p{u)Af, + p'{u)\V/h"{u) and expanding the square, we arrive at 



C"(0) = - j (^{C^,Ka{u)p{u) + h'\u)p{uf) (A^^ 

(17) + (cRKa(u)^ + 2p(u)/i'(u)) Ae|Ver + ^|Ver)da; 

= - f ((Crk + l)h’'{u)p{ufifl + (Crk + 2)p{u)p’{u)fLfG + T'iufh"{u)~^f^)dx, 


where we have employed the identity a{u) = p{u)h''{u) and the abbreviations fa = |V^| 
and fL = Af,. 
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We apply now the method of systematic integration by parts [H] . The idea is to identify 
useful integration-by-parts formulas and to add them to G"(0) without changing the sign 
of G"(0). The hrst formula is given by 


(18) 


div (ri(ti)(vy - Ap) ■ V5)di = / 


'an 


where ri(M) < 0 is an arbitrary (smooth) scalar function which still needs to be chosen, 
and I is the unit matrix in The left-hand side can be expanded as 

h''{u) 


Ve' - Ap)Ve + ri(n)V^e : - AP)J dx 

\ghg - TTrt^^L^G + ~ dx, 


h"{u) 


h"{u) 


where we have set ^ghg = and The boundary integral in flT8|) 

becomes 


/an 


ri(w) xV(|Ven - Aeve )-uds = - ri(n)V(|Ven ■ uds > 0, 


Ian 


since ri(M) < 0, ■ 1 / = 0 on dO,, and it holds that V(|V^p) • z/ < 0 on for all smooth 

functions satisfying ■ u = 0 on dQ [THl Prop. 4.2], Here we need the convexity of fl. 
Thus, the hrst integration-by-parts formula becomes 


(19) 


Ji-.= 


h"{u 

The second formula reads as 


\ghg - + ^ 1 ( 11 )^^ - Ti{u)^l ] dx >0. 


h"{u) 


( 20 ) 


0 = / div(r2(n)|V^pV0da: 


n{u) 

h"{u) 


+ ‘^^ 2 {u)^ghg + ) dx =: J 2 , 


where r 2 is an arbitrary scalar function. The goal is to hud functions ri(M) < 0 and r 2 (M) 
such that G"(0) < G"(0) -|- Ji -|- J 2 < 0. 

According to [15], the computations simplify if we introduce the variables and 
satisfying 

{d - l)eG^S = iGHG - la^Gi + d{d - 1)^1 + a- 


The existence of follows from the inequality 


d 


(I = > hAJ )2 + 


d 


VCV=8V{ A{V Irf , j,, 
ve ~irj + 


d d — 1 

which is proven in [JS] Lemma 2.1]. Then 

(21) G"(0) < G"(0) -|- Ji -|- J 2 = — f {ai^l + a2^L^G + + (^6^s)dx, 
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where 


( 22 ) 


^ (Crk + l)h"(u)^(uy + ^ j ri(u), 

02 = (Crk + 2)/4uj/y(uj +( 1 - 2 ) - 

4{nr-r'4u) , _ ,A£iM 

h''{u) 


as — 


h''{u) 


a 4 = -(h-l) 


a5 = -Ti{u), ae =-d{d-l)Ti{u). 


0 + 1) r 2 («), 

+ 2 r 2 (M)^ , 


The aim now is to determine conditions on oi,... ,06 such that the polynomial P{^) = 
+ ®4^s^G + IS nounegative as this implies that G"( 0 ) < 0 . 

In the general case, this leads to nonlinear ordinary differential equations for Ti and r 2 
which cannot be easily solved. A possible solution is to require that the coefficients of the 
mixed terms vanish, i.e. 02 = 04 = 0 , and that the remaining coefficients are nonnegative. 
The case d = 1 being simpler than the general case (since Ji is not necessary), we assume 
that d > 1. Then 04 = 0 implies that T[{u)/h''{u) = —2T2{u). Replacing T[{u)/h"{u) by 
— 2 r 2 (n) in 02 = 0 gives 

r 2 (M) = 

On the other hand, replacing r 2 (M) by —T[{u)/{2h''{u)) in 02 = 0, we hnd that 

T[{u) = -^(Ork + 2)n{u)n'{u)h''{u) 

or, after integration, 

ri(M) = -^(Ork + 2) f n{v)iJ,'{v)h"{v)dv. 

JuQ 

These functions have to satisfy the conditions 

oi > 0 or ^ V i(m) > -(Crk + l)h''{u)i4{uf, 

03 > 0 or (Crk + 2)^{u)^"{u) + (Crk - l)h'(“)^ < 0 , 

05 > 0 or ri(M) < 0 for all u, 

Note that Oi > 0 and 05 > 0 correspond to ffTSj) and (HI, respectively. This shows that 

F(0 > 0 for all ^ G and ^"(O) < 0. 

If G"(0) = 0, the nonnegative polynomial P, which depends on x G O via has to 
vanish. In particular, 03 ^^ = oslVul'^ = 0 in O. As 03 > 0 by assumption, u{x) = const, 
for X E fl. This contradicts the hypothesis that u is not a steady state. Consequently, 

G"(0) < 0, and we hnish the proof by setting b{u) = —ri(o). □ 
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4. Porous-medium equation 

The results of the previous section can be applied in principle to the Runge-Kutta scheme 
for the porous-medium or fast-diffusion equation 

(23) dtu = A(u^) in n, t > 0, ■ u = 0 on 912, n(0) = 

where f3 > 0. It can be seen that conditions (ll4l) - ffT6l) are not optimal for particular 
entropies. This is not surprising since we have neglected the mixed terms in the polynomial 
in fl 2 T]) (i.e. 02 = 04 = 0 ) which is not optimal. In this section, we make a different approach 
by making an ansatz for the functions Ti and r 2 , considering both zeroth-order and hrst- 
order entropies. 

4.1. Zeroth-order entropies. We prove the following result. 

Theorem 8. Let 12 C 6e convex with smooth boundary. Let {u^) be a sequence of 
(smooth) solutions to the Runge-Kutta scheme d^D for fl2^ . Let the entropy be given by 
H[u] = Q!“^(a!-|-1)“^ with a > 0, let k E N, and let be not the constant steady 

state of (l2^ . There exists a nonempty region Rnid) C (0, 00 )^ and > 0 such that for 
all (a,/3) ERo{d) and 0 < r < 

H[u^] +tP [ < H[u^-\ keN. 

In one space dimension, we have 

implicit Euler: Ro(l) = (0, cxo)^, 

Runge-Kutta of order p >2 : Rq{I) = {(a, (3) E (0, 00 )^ : —2 < a — (I < 1}, 

explicit Euler: Ro(l) = {(«)/9) E (0, 00 )^ : —1 < a —/3 < 1}. 

For the implicit Euler scheme, the theorem shows that any positive values for {a, ft) 
is admissible which corresponds to the continuous situation. For the Runge-Kutta case 
with Crk = 1; our condition is more restrictive. As expected, the explicit Euler scheme 
requires the most restrictive condition. The set Ro{d) is illustrated in Figured] for d = 2 
and d = 10 . 

Proof. Since A; G N is hxed, we set u := u^. We choose the functions 

Fi(n) = F 2 (n) = C 2 / 3 V^- 2 “-b 

It holds h"{u) = and /!(«) = ftu^~^. Then the coefficients in (|22|) are as follows: 

Oi = /3 ^((C'rk -M) (1 - 

02 = /3^((C'rk + 2){P - a) (1 - i)(2/3 -a - l)ci - (| + l)c2)n^^“^“"\ 

03 = /9"((/9 - - {2(3 -2a- l)c2)n2^-"“-^ 

04 = —ft‘^{d — l)((2/5 — a — l)ci -|- 2 c 2 )n^^“^““^, 

Os = —06 = —(3‘^d{d — . 
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Figure 1. Set Ro{d) of all {a, (3) for which the zeroth-order entropy is 
dissipating. Left column; d = 2, right column; d = 10. Top row; explicit 
Euler scheme with Crk = 2, middle row; implicit Euler scheme with Crk = 
1, bottom row; Runge-Kutta scheme of order p >2 with Crk = 0. 


Introducing the variables pj = for j G {G, L, R, S'}, we can write fl^ as 


^2/3+a 


G"(0) < G"(0) + + J 2 = -/32 
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where Q{r]) = birjl 621 /LI/G + hVc + biVsVG + ^Vr + 
with coefficients 


bi — (Crk + 1) + (1 ~ 

62 = (C*RK + 2)(/3 — «) + (! — ^)(2/3 — a — l)ci ~ (f + l)c 2 , 
h = W - a)^ - (2/3 -2a- l)c 2 , 

64 = -{d - l)((2/9 - a - l)ci + 2 C 2 ), 

65 = —Cl, = —d{d — l)ci. 


We need to determine all (a,/3) such that there exist Ci < 0, C 2 G M such that Qijf) > 0 
for all 7] = {rjc, ? 7 l, i/r, rjs). Without loss of generality, we exclude the cases 61 = 62 = 0 and 
&4 = fee = 0 since they lead to parameters {a, (3) included in the region calculated below. 
Thus, let fei > 0 and fee > 0. These inequalities give the bound — (C'rk + 1)/(1 —1/d) < ci < 
0. Thus, we may introduce the parameter A G (0,1) by setting ci = —A(Crk + 1)/(1 —1/d). 
The polynomial Qir]) can be rewritten as 


Qiv) = fei 

> Vg 


V 4fee AbJ 




+ 1 ) 


dfeifee 


2 

+ hvl + Vg 


R{c2;X,a,/3), 



4fee 


where R{c 2 ] X,a, (3) is a quadratic polynomial in C 2 with the nonpositive leading term 
—d^(4 —3A) + 4(2 —3A)d —4. The polynomial i?(c 2 ; X,a,{3) is nonnegative for some C 2 if and 
only if its discriminant 4d^A(l — A)S'(A; a, (3) is nonnegative. Here, S'(A; a, (3) is a quadratic 
polynomial in A. In order to derive the conditions on {a, (3) such that S{X]a,(3) > 0 for 
some A G (0,1), we employ the computer-algebra system Mathematica. The result of the 
command 

Resolve [Exists[LAMBDA, S[LAMBDA] >= 0 && LAMBDA > 0 
&& LAMBDA < 1], Reals] 

gives all {a, (3) G such that there exist ci < 0, C 2 G M such that Qir]) > 0. The interior 
of this region equals the set Ro{d), defined in the statement of the theorem. This shows 
that G"(0) < 0 for all {a, (3) G i?o(d). 

If G"(0) = 0, the nonnegative polynomial Q has to vanish. In particular, birjl = 0. 
li rjL = 0 in H, the boundary conditions imply that u is constant, which contradicts our 
assumption that u is not the steady state. Thus fei = 0. Similarly, fe 2 = fes = fe 4 = 0. 
This gives a system of four inhomogeneous linear equations for (ci, C 2 ) which is unsolvable. 
Consequently, G"(0) < 0. 

The set Ro{d) is nonempty since, e.g., (1,1) G Ro{d). Indeed, choosing ci = —1 and 
C 2 = 0, we find that Q{r]) = (Grk + ^)hi + hi + d(d - l)r]g > 0. 
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In one space dimension, the situation simplifies since the Laplacian coincides with the 
Hessian and thus, the integration-by-parts formula ffT^ is not needed. Then (see fl2U]) I 

G"(0) = G"(0) + Ji = -(5^ f 

Jn 

where 

oi = Crk + 1, 02 = (C'rk + 2)(/3 — a) — 3 c 2 , 03 = (/3 — a)^ — (2/3 — 2a — l)c 2 . 

The polynomial P{^) = + ® 2 i/ + 03 ) with y = ^l/^g nonnegative if and only if 

Oi > 0 and 40303 — O 2 > 0 , which is equivalent to 

(24) - 9c^ + 2((C'rk -2){a-(3) + 2(C'rk + l))c 2 - Cl^{a - > 0. 

This inequality has a solution C 2 G M if and only if the quadratic polynomial has real roots, 
i.e. if its discriminant is nonnegative, 

0 < ((C'rk - 2)(a - /3) + 2 (Crk + 1))' - 9Cl^{a - 
= 4(Crk + 1) (—(2Crk — l)(a — /3)^ -|- (Crk — 2){a — (3) + (Crk + 1)) • 

The polynomial — (2Crk — 1)^^ + (Crk — 2)z -|- (Crk + 1) with z = a — (3 is always 
nonnegative if Crk = 0 (implicit Euler). For Crk = 1 and Crk = 2 , this property holds if 
and only if —(Crk + 1)/(2Crk — 3) < a — (3 < 1. This concludes the proof. □ 


4.2. First-order entropies. We consider the one-dimensional case and first-order en¬ 
tropies with f{u) = u°‘P, a > 0 . 

Theorem 9. Let Q <Z M. be a bounded interval. Let {u^) be a sequence of (smooth) solutions 
to the Runge-Kutta scheme (EJ of order p > 2 for fl2^ in one space dimension. Let the 
entropy be given by F[u] = J^(u°‘^‘^)ldx with a > 0, let k G be fixed, and let be not 
the constant steady state of (12^ . There exists a nonempty region Ri G [0, 00 )^ and > 0 
such that for all (a, (3) G Ri, there is a constant Ca,g > 0 such that for all 0 < r < , 

F[u^] + rC„,^ [ < F[u^-^], keN. 

Jn 

Figure |2] illustrates the set Ri. The set of admissible values {a, (3) for the continuous 
equation is given by {—2 < a — 2/3 < 1} (the borders of this set are depicted in the figure 
by the dashed lines). 


Proof. First, we compute C'(0) according to Theorem [21 

C'(0) = -a [ 

Jn 

We show that C'(0) is nonpositive in a certain range of values {a, (3). We formulate C'(0) 
as 


C'(0) = 




4 
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Figure 2. Set of all {a, (3) for which the discrete hrst-order entropy for 
solutions to the one-dimensional porous-medium equation is dissipating. The 
continuous hrst-order entropy is dissipated for —2 < a —2/3 < 1. The borders 
of this set is indicated in the hgure by dashed lines. 

where = u^/u, ^2 = u^x/u. We employ the integration-by-parts formula 

0 = [ {u°‘'^^~^ul)xdx = [ {{a + (3 - ?,^l^2)dx =: J. 

Jo. Jn 

Therefore, 

G'(0) = G'(0) - ^cJ = -31 j 

where 

-P(0 = ((ci — 2)(/3 — 1) -|- (a -1- /3 — 4)c).^^ + {a + 2(3 — A + 3c).^^.^2 + ‘^^ 2 - 
This polynomial is nonnegative if and only if 

8((a - 2)(/3 - 1) + (a + /3 - 4)c) - (a + 2/3 - 4 + ?,cf > 0, 
which is equivalent to 

gi^c) := -9c^ + 2(a - 2/3 - 4)c - (a - 2/3)^ > 0. 

The maximizing value c* = {a — 2(3 — 4)/9, obtained from g\c) = 0, yields 

g{c*) = -^(a - 2/3 - l)(a - 2/3 + 2) > 0 

and consequently G'(0) < 0 if —2 < a — 2(3 < 1. This condition is the same as in [6l 
Theorem 13] for the continuous equation. 

Next, we turn to the proof of G"(0) < 0. The proof of Theorem [3] shows that 

G"(0) = (f - (f - 
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^{u^)xx)^^dx. 

We integrate by parts in the last term and use {(3u^~^{u^)xx)x = 0 on dil: 



X (®l'Cl + 0 . 2 ^ 1^2 + OS'ClCs + 0,4^i^2 + 0.5^1^2^3 + 0 - 6^2 + 0 . 7 ^ 3 ) dx, 
where =Ux/u, ^2 = Uxx/u, ^3 = Uxxx/u, and 

Oi = (/3 — 1) (2C*Rj<;cr^/3 — 3 C*rrci^ + 2(^(3“^ — 2(5C*p{r + 3)q!/3 + (15C*rx + 4)ci 
+ 2/3^ - 14/32 + 4(3Crk + 7)(3 - 2(9C'rk + 8 )), 

02 = (/3 — l)(4C'RKa^ + (8 Crk + 7)a(3 — (32Crk + 9)a + 12f3‘^ — 2(8Crk + 25)/3 
+ 6(8Crk + 7)), 

03 = C'rkci^ + 2a(3 — (5Crk + 2 )a + 4(C'rk + l)/?^ — 2(5 C'rk + 8)/3 + 12(C'rk + 1), 

04 = 2(/3 — 1)(2(4Crk + l)a + 9/3 — (16 Crk + 13)), 

05 = 2(2Crk + l)a + 4(2 C'rk + 3)/3 — 16(Crk + 1), 

06 = 2 — a, 07 = 2 (Crk + 1 ). 

We employ three integration-by-parts formulas: 

0 = / (u'^+^^-^ul^Ux) Jx = [ m“+2/3-22/3-5)^2^2^ 2 ^ 1^26 + ^1)dx =: Ji, 

Jn Jn 

0 = [ {u°'+'^^~^Uxxul)^dx = [ u^+'^^~^{{a + 2l3-Q)^f ^2 + ^U3 + ^^l^i)dx =: J 2 , 

Jn Jn 

0 = f Jx = [ m"+2/3-2^^^ ^2/3- 7)^4 + =: J 3 . 

Jn Jn 

Then 

G"(0) = G"(0) - la^f3\ciJi + C 2 J 2 + C 3 J 3 ) = -^«'/3' [ u^+^f^-^P{0dx, 

8 ^ Jn 

where -P (0 = + ^5^1^2^3 + ^ 6^2 + ^7^3; 

and the coefficients are given by 

61 = cti -|- (cr -|- 2/3 — 7 )c 3 , 62 = 02 + (cr + 2/3 — 6 )c 2 -|- 5 c 3 , 

^3 = 03 -|- C 2 , 64 = (34 -f (cr -|- 2/3 — 5)ci -|- 3 c 2 , 

^5 = 05 -f- 2 ci, fcg = Oe + Cl, 

67 = (Xj. 

Choosing ci = —Og, we eliminate the cubic term Furthermore, setting, x = ^ 2/^1 and 
y = we can write the polynomial P as a quadratic polynomial in {x, y): 

Q{x, y) = ^iP(0 = hi + b 2 X + bsy + b^x'^ + b^xy + bjy"^. 
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The following lemma is a consequence of the proof of Lemma 2.2 in [16]. 

Lemma 10. The polynomial p{x, y) = A + Bx + Cy + Dx^ + Exy + Fy"^ with F > 0 is 
nonnegative for all {x, y) G if and only if 

(i) ADF -F^ >0 and A{ADF - F^) - B‘^F - C^D + BCF > 0, or 

(ii) ADF — F^ = Q and 2BF — CF = 0 and AAF — > 0. 

Note that in case ADF — F"^ = 0 and F ^ 0, we may replace 2BF — CL' = 0 by the 
condition 2BEF = CL^ = ACDF or (since L > 0) BE = 2CD. 

The hrst inequality in case (i), 

0 < 46467 — 65 = —(Crk + 1)(2Crk + + (2Crk + 2)(4Crk — 3)q!/3 + (9Crk + 9)a 

— 2Crk(4Crk + 3)/3^ + (8Crk + 12)/3 + (3Crk + 3)c2 — (12 Crk + 14), 

is linear in C 2 and provides a lower bound for C 2 : 

^ q77 ^-TTi 4“ 1)(2 C'rk + l)a^ — (2Cpi,K + 2)(4Crk — 3)q;/3 — (9Crk + 9)a 

oIgrk + tj V 

+ 2Crk(4Crk + 3)/3^ — ( 8 Crk + 12)/3 + (12Crk + 14) j =: c^. 

The second inequality in case (i) becomes 

0 < 61(46467 - 65 ) - 6267 - 6364 626365 = -50(Crk + l)c 3 +pi(a,/3,C2)c3 +P2(a,/4,C2), 

where pi and p 2 are some polynomials in a, /3, and C 2 . This quadratic expression in C 3 is 
nonnegative if and only if its discriminant is nonnegative, 

0 < -200(Crk + I)p2(a,/5,C2) - pi(a,/5, 02 )^ 

= -8(46467 - 6^) ( 25 c 2 +P3(a,/4)c2 +P4(a,/4)), 

where p 3 (a,/ 3 ) and p 4 (a,/ 9 ) are some polynomials in a and (3. The factor 46467 — h\ is 
positive, so we have to ensure that Ra, 0 {,C 2 ) = 25cl + ps^a, ( 3 )c 2 + P 4 {cx, (3) < 0 for some 
C 2 > C 2 . Therefore we must ensure that the rightmost root of Ra, 0 {c 2 ) is larger or equal 
than the lower bound for C 2 , i.e., —^ 3 ( 0 , (3) + ^/p 3 {ct, (3) — 100 p 4 (a,/ 3 ) > SOc^. For Crk = 1, 
the values (a, (3) for which there exists C 2 > C 2 such that Ra, 0 {c 2 ) < 0 is depicted in Figure 
121 In case (ii), we may immediately calculate C 2 and C 3 but this results in a region which 
is already contained in the first one. This shows that G"(0) < 0. 

If G"(0) = 0, the polynomial Q vanishes. Thus, either u^/u = .^1 = 0 or P{^) = 0 in hi. 
The hrst case is impossible since u is not constant in hi. As bj = qt = 2(Grk + 1 ) > 0 , the 
second case P{ff) = 0 implies that ^3 = 0. Hence, m is a quadratic polynomial. In view of 
the boundary conditions, u must be constant, but this contradicts our assumption. Hence, 
G"(0) <0. □ 
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5. Linear diffusion system 


We consider the following linear diffusion system: 

(25) - PiAmi = /i(M2 - Ml), dtU2 - P2^U2 = - U2), 

with initial and homogeneous Neumann boundary conditions, pi, p 2 , p > 0, and the 
entropy 

(26) H[u] = f h{u)dx = 

Jn 

where u = {ui,U 2 )- If the initial data is nonnegative, the maximum principle shows that 
the solutions to fl 2 ^ are nonnegative too. 


/ y^^Ui{\ogUi - l)dx, 
Jn -1 


Theorem 11. Let {u^) he a sequence of (smooth) nonnegative solutions to the Runge-Kutta 
scheme (Ej) for with Crk = 1 and p := pi = P 2 - Let the entropy H be given by fl26l) . 
Let k eN be fixed and let be not the steady state of ([2]) . Then there exists > 0 such 
that for all 0 < T < 


E 

2=1 


|Vm: 


k\2 




+ p(logM^ - logM^)(M^ - M^) ]dx < 


Note that we need equal diffusivities pi = p 2 and higher-order schemes (Crk = !)• These 
conditions are in accordance of [IB], where the continuous equation was studied. In order 
to highlight the step where these conditions are needed, the following proof is slightly more 
general than actually needed. 

Proof. We fix k E N and set u := u^. Let A[u] = {Ai[u], A2[u]) = (piAwi -|- p{u2 — 
Mi),P 2 Am 2 + p(mi — ^2)). Since A is linear, DA[u\{h) = A[h]. Thus, 

G"( 0 ) = - [ {CKKh'{ufA[A[u]] + A[u]'^h'\u)A[u])dx = -Gi - G2. 

Jn 

In the following, we set dih = dh/dui for i = 1, 2. We integrate by parts twice, using the 
boundary conditions Vwj • z/ = 0 and • z/ = 0 on 912, and collect the terms: 

Cl = Crk j {dih{u){piJ^Ai[u] + p{A2[u] - Ai[m])) 

+ d2h{u)[p2JAA2[u] + p{Ai[u] - A2[u]))^dx 

= Crk j {j3iMih{u)Ai[u]+ p2M2h{u)A2[u] 

+ p{dih{u) — d2h{u)){A2[u] — Ailuf^dx 

= Crk j (pi{dlh{u)Aui + dlh{u)\Vui\^){piAui+p{u2-ui)) 

+ P2{dih{u)Au2 + d2h{u)\'\/u2\‘^) (P 2 AM 2 + p(mi - U 2 )) 
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fi{d2h{u) - dih{u))[piAui - P2AU2 + 2 p{u 2 - ui)) jdx 


= C'rk J ypidih{u){Auiy + P2^2^(m)(Am2)^ + pidfh{u)Aui\Vuif 

+ P 2 ^ 2 ^(m)Am 2 | V'U 2 p + Pip{dlh{u){u 2 — Ml) + d 2 h{u) — dih{u))Aui 
+ P2p{djh{u){ui - U2) + dih{u) - d2h{u))Au2 + pipdlh{u){u2 - mi)|Vmi| 
+ P2pdlh{u){ui - M2)|VM2p + 2p^{d2h{u) - dih{u)){u2 - Ui))dx. 


Furthermore, 


G 2 = 


dlh{u)(piAui + p{u2 - ui)) + dlh{u)[p2Au2 + p{ui - U2)) jdx 
pldlh{u){Auiy + P 2 ^ 2 ^(m)(Am 2 )^ + 2pipd‘lh{u){u2 — Ui)Aui 


2p2pd2h{u){ui — U2 )Au2 + p^{dlh{u) + c ? 2 ^( m))(mi — U2Y jdx. 


Adding Gi and G 2 , we arrive at 
2 


<^"(0) = - ^ / (Pi (Crk + l)dih{u){Auiy + p‘fGnKdfh{u)Aui\Vui\‘^'^dx 

i=l 

- J (pip((Crk + 2 )dlh{u){u 2 - Mi) + Crk( 52 ^(m) - dih{u)))Aui 

+ P 2 p((Crk + 2)dlh{u){ui - U 2 ) + Crk(5i^(m) - d2h{u)))Au2 
+ pi/iCRKl9^h(M)(M2 - Mi)|VmiP + p2pG^^dlh{u){ui - M2)|VM2nc^a; 


- / p \y2{dih{u)-d2h{u)) + {dih{u) + d2h{u)){ui-u2)j{ui-u2)dx 

J 

= -h -h- /o. 

The idea of [H] is to show that each integral /*, involving only derivatives of order i, is 
nonnegative. In contrast to [18], we employ systematic integration by parts, which allows 
for a simpler and more general proof in our context. For the term J 2 , we use the following 
integration-by-parts formula: 

0= / div (m“^| = / ( — 2 m“^|-|-SmC^AmjI =: Jj. 

Then, for £ > 0, 


I 2 - c'^p^Ji - e ^\'Vui\‘^dx 


i=l 


i=l 
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(^(C'rk + ^{AUiY 


(3c + Crk)^^ Aui\VUi\‘^ + (2c — e)u^ '^\Vui\^'^dx. 


The integrand defines a quadratic polynomial in the variables Aui and iVwjp and is non¬ 
negative if its discriminant satisfies 4(2c — e)(CRK + 1) ~ (3c + Crk)^ > 0. It turns out 
that this inequality holds true for Crk G {0,1} if we choose c = 2/3 and e > 0 sufficiently 
small. When Crk = 2, we can show only that /2 > 0 which is not sufficient to prove that 
G"(0) < 0 (see below). We conclude that 

2 „ 

(27) uf\Vui\^dx. 

i=i 


Integrating by parts in Ii in order to obtain only first-order derivatives, we find after 
some rearrangements that 


Ji =/i / (ai|V logMil^-I-a 2 V logMi • V logM 2 + oalV logM 2 |^)da;, where 

Jn 

oi = 2 pi(Crk'Ui + M2), 0,3 = 2p2(C'rkM2 + mi), 

02 = —(C'rk(Pi + P 2 ) + 2 ^ 2 )^! — (Crk(Pi + P 2 ) + 2Pi)M2. 


The integrand is nonnegative if and only if doiOa — > 0 for all (mi,M 2 )- We compute: 

Crk = 0 : doios -al = -4(piM2 - P2Ml)^ 

Crk = 1 : 40103 - 02 = (pi - P2)(pi(m? + QU1U2 + 9 ul) - P 2 ( 9 m^ -h 601^2 + uD), 
Crk = 2 : 401O3 - 02 = - 4 (pi(mi 2^3) - P 2 ( 2 mi + U2 )). 


Thus, 4 oi 03 — O 2 > 0 is possible only if pi = p 2 and Crk = 1. 
Finally, we see immediately that the remaining term 



|^2(logMi - logM2)(Ml - M 2 ) + 



(Mi - M 2 ) 


dx 


is nonnegative. This shows that G"(0) < 0. If G"(0) = 0, we infer from (HT} that 
Mj = const., but this contradicts our hypothesis that Mj is not a steady state. □ 


6 . The Derrida-Lebowith-Speer-Spohn equation 

Consider the one-dimensional fourth-order equation 

(28) dtu =-{u(logu)^^)^^ mQ,t>o, m(0) = m° 

with periodic boundary conditions. This equation appears as a scaling limit of the so-called 
(time-discrete) Toom model, which describes interface fluctuations in a two-dimensional 
spin system [9] . The variable m is the limit of a random variable related to the deviation of 
the spin interface from a straight line. The multi-dimensional version of (|28|) models the 
eectron density m in a quantum semiconductor, und the equation is the zero-temperature, 
zero-field approximation of the quantum drift-diffusion model [13]. For existence results 
for (|28|) . we refer to [15] and references therein. 
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To simplify our calculations, we analyze only the logarithmic entropy H[u] = u(logu — 
l)dx. It is possible to verify condition (| 6 ]) also for entropies of the form f^u‘^dx, but it 
turns out that only sufficiently small a > 0 are admissible (about 0 < a < 0.15...) and 
the computations are very tedious. Therefore, we restrict ourselves to the case a = 0. 


Theorem 12. Let (u^) be a sequence of (smooth) solutions to the Runge-Kutta scheme 
(E]) with Crk = 1 for (1^ . Let the entropy be given by H[u] = J^u(\ogu — l)dx, let 
k E N be fixed, and let be not a steady state. Then there exists > 0 such that for all 
0 < r < 

+ Tq u{logu)ldx + T 'o{logu)l„,dx < q ^ 0.0045. 

Jq Jq, 


Proof. First, we observe that G'( 0 ) = — J^{u{logu)xx)xx logudx = — Jq u(logu)^,j.dx. With 
A[u] = {u{\ogu)xx)xx and DA[u]{h) = {hxx - 2{logu)xhx + {logu)lh)^^, we can write 
G"(0) = —Jq according to ([H]) as 


G"(0) = - / ( logu{A[u]xx - 2(logu)xA[u]x + (logu)lA[u ])+ -A[u]Adx 

.In ^ U / 


(\ogu)xx{A[u]xx - 2(\ogu)xA[u]x + (\ogu)lA[u]) + -A[uf)dx 

u 

1 


I ((oxxxx T 2(VxVxx')x T T '\dx, 

Jn ^ u / 

where we have integrated by parts several times and have set v = logn. Then A[u] = 

u(vlvxx + 2vxVxxx + pL + Oxxxx) and, with the abbreviations 6 = W, • • •, ^4 = Vxxxx, 


G"(0) 


“ j + 4'Ci'C2'C4 + §.^ 1.^3 + lO.^i.^^'Cs 

+ 8'Ci'C3'C4 + 3.^2 + ^^2^4: + 2ff^dx. 


We employ the following integration-by-parts formulas: 


0 = / ( uvl ) xdx = f u ( f ^ + 7 f%)dx =: Ji , 

J ft J ft 

0 = [ (uVxxvDxdx = [ u{f% + =■ J 2 , 

•Jn Jn 

0 = / ( uVxxxvt)xdx = / + ^3^4 -F 4 : f , ff , 2 f 3 )dx =: J3, 

Jn Jn 

0 = / (uvlxvDxdx = [ m(^i ^2 + + 3^i^2)c^^ =: -^4, 

Jn Jn 

0 = / (uVxxVxxxvDxdx = / + ^?6^4 + =: Jg, 

Jn Jn 

0 = / (uvlxxVx)xdx = / + 2^i^3^4 + f 2 Qdx =: Je, 

Jn Jn 
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0 = / iuv''^^V:,):,dx = I + 36'^26 + ^t)dx =: J 7 , 


0 / {uVxxxXxx)xdx / u(^^i^ 2^3 ^2^4:)dx . Jg. 


Then 


G"{0) — G"{0) — 4 CiJi — — u(ai^f + a 2 ^i ^2 + 0,3^l^3 + 0,4:^l^2 + ®5'Cl'C4 

i=i ^ 


+ + 07^1^2 + ®8^1^2^4 + <39^1^3 + ®10^1^2^3 + 011^1^3^4 + Ol2^2 

+ ai 3 'C 2 'C 4 + aM^ 2 ^i + 0 , 15 ^ 1 '^ dx, 


- 2 t 2 


where 


oi = 4ci, 

04 = 2 T 20c2 T 4c4, 
07 = 5 + 12 c 4 + 4 c 7 , 


0,2 — 28ci + 4 c2, 
05 = 4c3, 

08 = 4 + 4 c 5 , 


Olo — 10 + 8C5 + 12c7 + 4 c 8 , Oil — 8 T 8cg, 


fll3 — 5 + 4 c8, 


Oi4 — 4 c6 + 8 C 8 , 


03 = 4c2 + 4c3, 

Og = 8 + I6C3 + 8C4 + 4 c5, 

O9 = 8 + 4c5 + 4 cg, 

O12 = 3 + 4 c 7 , 

015 = 2 . 


Next, we eliminate all terms involving ^4 by formulating the following square: 

2 


G"(0) = - 


u 


In 


Oi5 ( .^4 + 77 -^^! 


2 a 


0-8 ^2c I C C I ®13 c2 

S1S2 + 77-a 43 + 7 ^-S2 


15 


2 o 


15 


2 o 


15 


2 a 


15 


Oi - 


4o 


15 


^ 1 ^ + a 2 


f 05as\ Q f ®5aii\^5c 


04 


09 


4o 


15 




C C < I 08^13^ ^2 <i 3 


O 


4o 


15 


11 \c'^c2 I / OllfllS \ * ^2c I / 

^ 1^3 + ( ®10 - ^ - 1 ^ 1 ^ 2^3 + ( 012 


4o 


15 


'^2 + 0146^3 


dx. 


We eliminate all terms involving ,^3 and set the corresponding coefficients to zero. From 
Oi 4 = 0 we conclude that Cg = —2c8. Furthermore, 


o 


09 

Oio - 
Og — 
03 - 


11 


40i5 

Q 11013 

20i5 

OsOii 

20i5 

05011 


2 a 


= 0 
= 0 
= 0 
= 0 


gives C 5 = 8 c 8 - 6 c 8 , 

20 2 8 

gives C7 = -yC8 + -C8, 

gives C 4 = — 2 c 3 — lOcg + 16cg — Scg, 

gives C 2 = C 3 - 4 c 3 C 8 . 


15 
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By these choices, we obtain 


bi 2 


■— 0,12 — 


4ai5 



This quadratic polynomial in cg admits its maximal value at Cg 
bi 2 = 20/129. The integral can now be written as 


17/172 with value 


where 


G"(0) < — / + ^7^1^2 + ^12^2)*^^) 


61 = ai - 


4a 


= 4ci - 2 C 3 , 


15 


b 2 = 02 - = 28ci - 32 C 3 C 8 + ScsCg, 


64 — 04 — 


2ai5 

Os 


4a 


15 


05^13 

2 a 


= 7c3 — 84c3C8 — 128c 8 + 128c8 — 40c8 + 4c8, 


15 


cir^is ^ a ^ a a r 448 9 70 

b-j = aj -= —24c 3 — 244co H-Co-Cg. 

^ ^ 2ai5 ® 3 ® 3 ^ 

If 64 = 262 & 12/&7 + &7/(46i 2), we can write the integal as the sum of two squares, noting 
that bi 2 is positive. 


G"(o) < - h^yijdx. 

The expression 6467 — 262&12 — &?/(45i2) = 0 dehnes a polynomial in (ci, C 3 ) which is linear 
in Cl- Solving it for ci gives 

449307 n 741681 , 35780649411 34135130165539 

^ “ 175 ^ 2150 ^ 2393160700 ^ 163091166664200‘ 

It remains to show that p{c^) := 61 — b\bi 2 /b‘^, which is a polynomial of fourth order in C 3 , 
is positive. Choosing Cg = —0.029, we hnd that ^(cg) ~ 0.0045 > 0. This shows that 


G"'(0) <-g(c;) [ u^fdx = -q{cl) [ u{logu)ldx < 0. 

Jn Jn 

Finally, if G"(0) = 0, we infer that u is constant which is excluded. Therefore, G"(0) < 0, 
which ends the proof. □ 


7. Numerical examples 

The aim of this section is to explore the numerical behavior of the second-order derivative 
of the function G(r), dehned in the introduction, for the porous-medium equation fl23|) in 
one space dimension. The equation is discretized by standard hnite differences, and we 
employ periodic boundary conditions. The discrete solution u/ approximates the solution 
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u{xi,t^) to fl^ with Xi = iAx, A = kr, and Ax, r are the space and time step sizes, 
respectively. We choose the Barenblatt prohle 


(29) 



^0 


max 



/3-1 

2/3(/3 + l) ) 


0 < X < 1, 


where 

to 0.01, C 2i3{fl+l) jW+i) ’ “ 4’ 

as the initial datum. Its support is contained in — xr, | + xr]] see Figure [3] (left). We 
choose the exponent (3 = 2. The semi-logarithmic plot of the discrete entropy Hd[u^] = 
Ylii=ohA )°'with a = 5 versus time is illustrated in Figure [3] (right), using the implicit 
Euler scheme with parameters r = 10“'^ and the number of grid points N = 1/Ax = 64. 
The decay is exponential for “large” times. The nonlinear discrete system is solved by 
Newton’s method with the tolerance tol = 10“^®. We have highlighted four time steps 
ti at which we will compute numerically the function G{t) for the following Runge-Kutta 
schemes: 


explicit Euler scheme: 


implicit Euler scheme: 


second-order trapezoidal rule: 


third-order Simpson rule: 



-rA[u% 

-^iA[u^] + A[u^-^]), 

~{A[u'^] + AA[{u^ + n"-')/2] + A[u^-^]). 
6 



Figure 3. Left: Evolution of the intial datum (1291) for 13 = 2 at various time 
steps U, i = 0,1,2, 3. Right: Semi-logarithmic plot of the discrete entropy 
Hd[u^] versus time. 

We set as before u := u^, v{t) := and compute G{t) = Hd[u] — Hd[v{T)] and 
the discrete second-order derivative d^G of G (using central differences). The result is 
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presented in Figure 01 As expected, the discrete derivative is negative on a (small) 
interval for all times tj, i = 1, 2, 3. We observe that d^G is even slightly decreasing, but 
we expect that it becomes positive for sufficiently large values of r. Clearly, the values 
for d^G tend to zero as we approach the steady state (see Remark 0)). This experiment 
indicates that from Theorem [T]is bounded from below bv r* = 3 ■ ID for instance. 




Figure 4. Numerical evaluation of the discrete version of G"{t) for various 
Runge-Kutta schemes at the time steps U. Top left: explicit Euler scheme; 
top right: implicit Euler scheme; bottom left: implicit trapezoidal rule; bot¬ 
tom right: Simpson rule. 


In order to understand the behavior of G{t) in a better way, it is convenient to study 
the discrete version of the quotient 


(30) 


Q(r) 


G"{r) 


Indeed, the analysis in Section 0] gives an estimate of the type G"(0) < —G ^u^dx 

for some constant C > 0. Thus, we expect that for sufficiently small r > 0, Qij) is bounded 

























































^ X 10 



'T X 10 




T X 1 0 


T X 1 0 


Figure 5. Numerical evaluation of the discrete version of Q{t), defined 
in flHU]) . for various Runge-Kutta schemes at the time steps ti. Top left: 
explicit Euler scheme; top right: implicit Euler scheme; bottom left: implicit 
trapezoidal rule; bottom right: Simpson rule. 


from above by some negative constant. This expectation is conhrmed in Figured In the 
examples, Q{t) is a decreasing function of r, and Q(0) is decreasing with increasing time. 

All these results indicate that the threshold parameter in Theorem [1] can be chosen 
independently of the time step k. 


References 

[1] D. Bakry and M. Emery. Diffusions hypercontractives. In: Seminaire de prohahilites XIX, 1983/84, 
Lect. Notes Math. 1123 (1985), 177-206. 

[2] M. Bessemoulin-Chatard. A finite volume scheme for convection-diffusion equations with nonlinear 
diffusion derived from the Scharfetter-Gummel scheme. Numer. Math. 121 (2012), 637-670. 

[3] S. Boscarino, F. Filbet, and G. Russo. High order semi-implicit schemes for time dependent partial 
differential equations. Preprint, 2015. https://hal.archives-ouvertes.fr/hal-00983924. 

[4] M. Bukal, E. Emmrich, and A. Jiingel. Entropy-stable and entropy-dissipative approximations of a 
fourth-order quantum diffusion equation. Numer. Math. 127 (2014), 365-396. 
















































ENTROPY-DISSIPATING SEMI-DISCRETE RUNGE-KUTTA SGHEMES 


29 


[5] C. Gances and C. Guichard. Numerical analysis of a robust entropy-diminishing finite-volume scheme 
for parabolic equations with gradient structure. Preprint, 2015. arXiv: 1503.05649. 

[6] G. Ghainais-Hillairet, A. Jiingel, and S. Schuchnigg. Entropy-dissipative discretization of nonlinear 
diffusion equations and discrete Beckner inequalities. To appear in ESAIM Math. Model. Numer. 
Anal, 2015. arXiv:1303.3791. 

[7] S. Ghristiansen, H. Munthe-Kaas, and B. Owren. Topics in structure-preserving discretization. Acta 
Numerica 20 (2011), 1-119. 

[8] K. Deimling. Nonlinear Functional Analysis. Springer, Berlin, 1985. 

[9] B. Derrida, J. Lebowitz, E. Speer, and H. Spohn. Fluctuations of a stationary nonequilibrium interface. 
Phys. Rev. Lett. 67 (1991), 165-168. 

[10] F. Filbet. An asymptotically stable scheme for diffusive coagulation-fragmentation models. Commun. 
Math. Sci. 6 (2008), 257-280. 

[11] A. Glitzky and K. Gartner. Energy estimates for continuous and discretized electro-reaction-diffusion 
systems. Nonlin. Anal. 70 (2009), 788-805. 

[12] E. Hairer, S.P. Nprsett, and G. Wanner. Solving Ordinary Dijferential Equations I. Springer, Berlin, 
1993. 

[13] A. Jiingel. Transport Equations for Semiuconductors. Lect. Notes Phys. 773, Springer, Berlin, 2009. 

[14] A. Jiingel and D. Matthes. An algorithmic construction of entropies in higher-order nonlinear PDEs. 
Nonlinearity 19 (2006), 633-659. 

[15] A. Jiingel and D. Matthes. The Derrida-Lebowitz-Speer-Spohn equation: existence, non-uniqueness, 
and decay rates of the solutions. SIAM J. Math. Anal. 39 (2008), 1996-2015. 

[16] A. Jiingel and J.-P. Milisic. A sixth-order nonlinear parabolic equation for quantum systems. SIAM 
J. Math. Anal. 41 (2009), 1472-1490. 

[17] A. Jiingel and J.-P. Milisic. Entropy dissipative one-leg multistep time approximations of nonlinear 
diffusive equations. To appear in Numer. Meth. Part. Diff. Eqs., 2015. arXiv: 1311.7540. 

[18] M. Liero and A. Mielke. Gradient structures and geodesic convexity for reaction-diffusion systems. 
Phil. Trans. Royal Soc. A. 371 (2013), 20120346 (28 pages), 2013. 

[19] H. Liu and H. Yu. Entropy/energy stable schemes for evolutionary dispersal models. J. Comput. Phys. 
256 (2014), 656-677. 

[20] E. Tadmor. Numerical viscosity of entropy stable schemes for systems of conservation laws 1. Math. 
Comp. 49 (1987), 91-103. 

Institute for Analysis and Scientific Gomputing, Vienna University of Technology, 
WiEDNER HAUPTSTRASSE 8-10, 1040 WiEN, AUSTRIA 
E-mail address: juengelStuwien. ac . at 

Institute for Analysis and Scientifig Gomputing, Vienna University of Technology, 
WiEDNER HAUPTSTRASSE 8-10, 1040 WiEN, AUSTRIA 
E-mail address: stefan.schuchnigg@tuwien.ac.at 


